TASR Analysis

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.3.6     ✔ purrr   0.3.4
## ✔ tibble  3.2.1     ✔ dplyr   1.1.2
## ✔ tidyr   1.2.0     ✔ stringr 1.4.1
## ✔ readr   2.1.2     ✔ forcats 0.5.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## 
## Attaching package: 'lubridate'
## 
## 
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
## 
## 
## corrplot 0.92 loaded
## 
## 
## Attaching package: 'ggpp'
## 
## 
## The following object is masked from 'package:ggplot2':
## 
##     annotate
## 
## 
## 
## Attaching package: 'flextable'
## 
## 
## The following object is masked from 'package:purrr':
## 
##     compose
## 
## 
## 
## Attaching package: 'maps'
## 
## 
## The following object is masked from 'package:purrr':
## 
##     map
## 
## 
## ℹ Google's Terms of Service: ]8;;https://mapsplatform.google.com<https://mapsplatform.google.com>]8;;
## ℹ Please cite ggmap if you use it! Use `citation("ggmap")` for details.
## 
## Attaching package: 'scales'
## 
## 
## The following object is masked from 'package:purrr':
## 
##     discard
## 
## 
## The following object is masked from 'package:readr':
## 
##     col_factor
## 
## 
## Loading required package: sp
## 
## Checking rgeos availability: TRUE
## Please note that 'maptools' will be retired by the end of 2023,
## plan transition at your earliest convenience;
## some functionality will be moved to 'sp'.
## 
## 
## Attaching package: 'raster'
## 
## 
## The following object is masked from 'package:dplyr':
## 
##     select
## 
## 
## rgeos version: 0.6-3, (SVN revision 696)
##  GEOS runtime version: 3.10.2-CAPI-1.16.0 
##  Please note that rgeos will be retired during October 2023,
## plan transition to sf or terra functions using GEOS at your earliest convenience.
## See https://r-spatial.org/r/2023/05/15/evolution4.html for details.
##  GEOS using OverlayNG
##  Linking to sp version: 1.6-0 
##  Polygon checking: TRUE 
## 
## 
## 
## Attaching package: 'rgeos'
## 
## 
## The following object is masked from 'package:dplyr':
## 
##     symdiff

Overview

This script reads in cleaned datafiles with stream temperature and invertebrate data produced by the invert_data_cleaning.Rmd and teton_data_cleaning.Rmd scripts and performs more involved analyses of stream temperature and invertebrate density/diversity over time. A few findings of interest from previous analyses are that:

* Average summer (July-Aug) stream temperatures have not changed significantly over time

* Average summer (July-Aug) air temperatures have increased over time

* Overall, invetebrate richness has not changed over time

* Invertebrate communities in snow-fed streams appear to be more variable than those in glacier-fed streams

To build on these findings, the main sections of this script are:

I. Temperature:

1. Overall trends in stream temperature over time - Time series analysis: Site-by-site and grouped by source

2. Trends in summer mean, max, min, and degree-days over time: Site-by-site and grouped by source

4. Modeling stream temperature using SNOTEL data: Site-by-site and grouped by source

II. Invertebrates:

1. Trends in richness, abundance, and density over time: Site-by-site and grouped by source

2. Modeling richness, abundance, and density using SNOTEL data: Site-by-site and grouped by source

III. Temperature and Invertebrates

1. Test for correlations between temperature metrics (mean/max/min summer; degree days): Site-by-site and grouped by source

Read in necessary files:

stream_temp <- read.csv("Temperature/cleaned_full_datasets/temps_hourly.csv") #hourly stream temp data
invert_density <- read.csv("invert_data/cleaned_csv/full_invert_densities.csv")%>%
  rename(site = Stream)#invert density data
invert_richness <- read.csv("invert_data/cleaned_csv/invert_richness.csv") #invert richness data
snotel <- read.csv("teton_snotel.csv") #snotel data
sources <- read.csv("source_info.csv")%>%
  rename(site = stream) #source info, rename for merge

### adding source info to stream temps and invert datasets: 

stream_source <- merge(stream_temp, sources, all = T)%>%
  mutate(date1 = ymd(date1)) #Add source info to hourly temperature data; re-code date as r date format
stream_source_daily <- stream_source%>%
  group_by(site, full_name, stream_code, date1, year, lat, long)%>%
  summarise(t_xbar = mean(temp_c, na.rm = T), 
            t_max = max(temp_c, na.rm = T), 
            t_min = min(temp_c, na.rm = T), 
            source = unique(source)) #calculate daily mean, min, and max temps for each site

## Warning: There were 10690 warnings in `summarise()`.
## The first warning was:
## ℹ In argument: `t_max = max(temp_c, na.rm = T)`.
## ℹ In group 41: `site = "cloudveil"`, `full_name = "Cloudveil"`, `stream_code =
##   "CV"`, `date1 = 2019-09-11`, `year = NA`, `lat = 43.72551`, `long =
##   -110.7961`.
## Caused by warning in `max()`:
## ! no non-missing arguments to max; returning -Inf
## ℹ Run ]8;;ide:run:dplyr::last_dplyr_warnings()dplyr::last_dplyr_warnings()]8;; to see the 10689 remaining warnings.

## `summarise()` has grouped output by 'site', 'full_name', 'stream_code',
## 'date1', 'year', 'lat'. You can override using the `.groups` argument.

density_source <- merge(invert_density, sources, all = T)
richness_source <- merge(invert_richness, sources, all = T)

I. Temperature:

1. Trends in summer mean, max, min, and degree-days over time

A. Site-by-site

First, calculate summer (july-aug) mean daily temp, avg daily max and min; plus degree days over 15c? for each site:

stream_summer <- stream_source_daily %>%
  mutate(mo = month(date1))%>% #add month column
  filter(mo == 7 | mo == 8)%>% #extract July and August
  group_by(site, full_name, stream_code, year, lat, long)%>% #group by site and year
  summarise(summer_xbar = mean(t_xbar, na.rm = T), 
            max_xbar = mean(t_max, na.rm = T), 
            min_xbar = mean(t_min, na.rm= T), #calculate mean temperature, mean daily max and min
            dd_2.5 = length(which(t_max >= 2.5)), 
            dd_5 = length(which(t_max >= 5)),
            dd_10 = length(which(t_max >= 10)), #calculate no. days with max temp >2.5, 5, and 10 deg
            source = unique(source)) #retain source col

## `summarise()` has grouped output by 'site', 'full_name', 'stream_code', 'year',
## 'lat'. You can override using the `.groups` argument.

Next, extract sites with at least 3 years of data:

under3 <- stream_summer %>%
  group_by(site)%>% #group by site
  filter(!is.na(year))%>% #remove rows with no year
  summarise(n_yrs = length(unique(year)))%>% #count num unique years for site
  filter(n_yrs < 3) #extract sites with <3yrs data

summer_3plus <- stream_summer %>%
  filter(!site%in%under3$site)

Trends in mean temp:

Nest data by site:

summer_nested <- summer_3plus %>%
  filter(!is.na(year))%>%
  nest(data = -site)

lm for each site:

summer.mean.lms <- summer_nested %>%
  mutate(model = purrr::map(data, ~lm(summer_xbar~year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
                     tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
  unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
  filter(term == "year")  #remove intercept term to just get info on year for each site
summer.mean.lms

## # A tibble: 11 × 7
## # Groups:   site [11]
##    site       data                term  estimate std.error statistic p.value
##    <chr>      <list>              <chr>    <dbl>     <dbl>     <dbl>   <dbl>
##  1 cloudveil  <gropd_df [3 × 12]> year    0.458     0.309      1.48   0.378 
##  2 delta      <gropd_df [4 × 12]> year    0.0246    0.0193     1.28   0.330 
##  3 grizzly    <gropd_df [4 × 12]> year    1.21      0.716      1.68   0.234 
##  4 mid_teton  <gropd_df [6 × 12]> year    0.0296    0.0280     1.05   0.351 
##  5 n_teton    <gropd_df [6 × 12]> year    0.255     0.460      0.554  0.609 
##  6 paintbrush <gropd_df [5 × 12]> year    1.06      0.510      2.08   0.129 
##  7 s_ak_basin <gropd_df [4 × 12]> year    0.0539    0.0159     3.40   0.0768
##  8 s_cascade  <gropd_df [3 × 12]> year   -0.0705    0.0549    -1.28   0.421 
##  9 s_teton    <gropd_df [6 × 12]> year    0.987     0.428      2.30   0.0826
## 10 skillet    <gropd_df [3 × 12]> year    0.370     0.0617     5.99   0.105 
## 11 windcave   <gropd_df [3 × 12]> year    0.119     0.0672     1.77   0.328

Nothing significant at 0.05 level; weakly significant increase in mean summer temp at S AK basin (p = 0.08)

Visualization:

summer_mean_withps <- merge(summer_3plus, summer.mean.lms)

pres.pal <- c("#C26ED6", "#F89225", "#76D96F")

summer.mean <- ggplot(summer_mean_withps, aes(x = year, y = summer_xbar, color = source))+
  geom_point()+
  geom_line()+
  geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
  theme_classic()+
  facet_wrap(~full_name, scales = "free")+
  scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
  geom_text_npc(aes(npcx = 0.9, npcy = 0.1, label = paste( "P = ",round(p.value, 2), sep = "")))+
  labs(x = "Date", y = "Mean daily July-Aug. stream temp., C", color = "Stream Source")+
  theme(legend.position = "bottom", 
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1), 
        axis.title = element_text(size = 11, face = "bold"), 
        axis.text = element_text(size = 8), 
        legend.text = element_text(size = 10), 
        legend.title = element_text(size = 11, face = "bold"), 
        strip.text = element_text(size = 10, face = "bold"))
summer.mean

## `geom_smooth()` using formula 'y ~ x'

## Warning: Removed 12 rows containing non-finite values (stat_smooth).

## Warning: Removed 12 rows containing missing values (geom_point).

## Warning: Removed 4 row(s) containing missing values (geom_path).

Mean daily max temp

summer.max.lms <- summer_nested %>%
  mutate(model = purrr::map(data, ~lm(max_xbar~year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
                     tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
  unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
  filter(term == "year")  #remove intercept term to just get info on year for each site
summer.max.lms

## # A tibble: 11 × 7
## # Groups:   site [11]
##    site       data                term  estimate std.error statistic p.value
##    <chr>      <list>              <chr>    <dbl>     <dbl>     <dbl>   <dbl>
##  1 cloudveil  <gropd_df [3 × 12]> year    0.532     0.328      1.62   0.352 
##  2 delta      <gropd_df [4 × 12]> year    0.0654    0.0816     0.801  0.507 
##  3 grizzly    <gropd_df [4 × 12]> year    1.66      1.10       1.51   0.271 
##  4 mid_teton  <gropd_df [6 × 12]> year   -0.0257    0.0485    -0.529  0.625 
##  5 n_teton    <gropd_df [6 × 12]> year    0.228     0.603      0.378  0.724 
##  6 paintbrush <gropd_df [5 × 12]> year    1.74      0.853      2.04   0.135 
##  7 s_ak_basin <gropd_df [4 × 12]> year    0.0774    0.108      0.713  0.550 
##  8 s_cascade  <gropd_df [3 × 12]> year   -0.130     0.133     -0.981  0.506 
##  9 s_teton    <gropd_df [6 × 12]> year    1.37      0.640      2.13   0.0999
## 10 skillet    <gropd_df [3 × 12]> year    0.365     0.0254    14.4    0.0442
## 11 windcave   <gropd_df [3 × 12]> year    0.158     0.101      1.56   0.364

Significant increase at Skillet, non-significant increases at all but S cascade and mid teton

Visualization:

summer_max_withps <- merge(summer_3plus, summer.max.lms)

summer.max <- ggplot(summer_max_withps, aes(x = year, y = max_xbar, color = source))+
  geom_point()+
  geom_line()+
  geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
  theme_classic()+
  facet_wrap(~full_name, scales = "free")+
  scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
  geom_text_npc(aes(npcx = 0.9, npcy = 0.1, label = paste( "P = ",round(p.value, 2), sep = "")))+
  labs(x = "Date", y = "Mean daily max. July-Aug. stream temp., C", color = "Stream Source")+
  theme(legend.position = "bottom", 
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1), 
        axis.title = element_text(size = 11, face = "bold"), 
        axis.text = element_text(size = 8), 
        legend.text = element_text(size = 10), 
        legend.title = element_text(size = 11, face = "bold"), 
        strip.text = element_text(size = 10, face = "bold"))
summer.max

## `geom_smooth()` using formula 'y ~ x'

## Warning: Removed 12 rows containing non-finite values (stat_smooth).

## Warning: Removed 12 rows containing missing values (geom_point).

## Warning: Removed 4 row(s) containing missing values (geom_path).

Mean daily min temp

summer.min.lms <- summer_nested %>%
  mutate(model = purrr::map(data, ~lm(min_xbar~year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
                     tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
  unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
  filter(term == "year")  #remove intercept term to just get info on year for each site
summer.min.lms

## # A tibble: 11 × 7
## # Groups:   site [11]
##    site       data                term  estimate std.error statistic p.value
##    <chr>      <list>              <chr>    <dbl>     <dbl>     <dbl>   <dbl>
##  1 cloudveil  <gropd_df [3 × 12]> year    0.400     0.273      1.46   0.382 
##  2 delta      <gropd_df [4 × 12]> year    0.0156    0.0152     1.02   0.413 
##  3 grizzly    <gropd_df [4 × 12]> year    0.849     0.493      1.72   0.227 
##  4 mid_teton  <gropd_df [6 × 12]> year    0.0639    0.0261     2.45   0.0708
##  5 n_teton    <gropd_df [6 × 12]> year    0.254     0.334      0.761  0.489 
##  6 paintbrush <gropd_df [5 × 12]> year    0.662     0.325      2.04   0.135 
##  7 s_ak_basin <gropd_df [4 × 12]> year    0.0398    0.0440     0.906  0.461 
##  8 s_cascade  <gropd_df [3 × 12]> year   -0.0249    0.0117    -2.14   0.279 
##  9 s_teton    <gropd_df [6 × 12]> year    0.667     0.292      2.28   0.0846
## 10 skillet    <gropd_df [3 × 12]> year    0.347     0.0963     3.61   0.172 
## 11 windcave   <gropd_df [3 × 12]> year    0.0993    0.0354     2.80   0.218

Non-significant increases at all sites except south cascade

Visualization:

summer_min_withps <- merge(summer_3plus, summer.min.lms)

summer.min <- ggplot(summer_min_withps, aes(x = year, y = min_xbar, color = source))+
  geom_point()+
  geom_line()+
  geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
  theme_classic()+
  facet_wrap(~full_name, scales = "free")+
  scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
  geom_text_npc(aes(npcx = 0.9, npcy = 0.1, label = paste( "P = ",round(p.value, 2), sep = "")))+
  labs(x = "Date", y = "Mean daily min. July-Aug. stream temp., C", color = "Stream Source")+
  theme(legend.position = "bottom", 
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1), 
        axis.title = element_text(size = 11, face = "bold"), 
        axis.text = element_text(size = 8), 
        legend.text = element_text(size = 10), 
        legend.title = element_text(size = 11, face = "bold"), 
        strip.text = element_text(size = 10, face = "bold"))
summer.min

## `geom_smooth()` using formula 'y ~ x'

## Warning: Removed 12 rows containing non-finite values (stat_smooth).

## Warning: Removed 12 rows containing missing values (geom_point).

## Warning: Removed 4 row(s) containing missing values (geom_path).

Degree Days

Days over 2.5 C:

summer.2.5.lms <- summer_nested %>%
  mutate(model = purrr::map(data, ~lm(dd_2.5~year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
                     tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
  unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
  filter(term == "year")  #remove intercept term to just get info on year for each site

## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `model = purrr::map(data, ~lm(dd_2.5 ~ year, data = .) %>%
##   tidy)`.
## ℹ In group 8: `site = "s_cascade"`.
## Caused by warning in `summary.lm()`:
## ! essentially perfect fit: summary may be unreliable

summer.2.5.lms

## # A tibble: 11 × 7
## # Groups:   site [11]
##    site       data                term   estimate std.error statistic p.value
##    <chr>      <list>              <chr>     <dbl>     <dbl>     <dbl>   <dbl>
##  1 cloudveil  <gropd_df [3 × 12]> year   6.00e+ 0  6.35e+ 0     0.945  0.518 
##  2 delta      <gropd_df [4 × 12]> year   4.43e+ 0  2.90e+ 0     1.53   0.266 
##  3 grizzly    <gropd_df [4 × 12]> year   6.1 e+ 0  5.16e+ 0     1.18   0.359 
##  4 mid_teton  <gropd_df [6 × 12]> year   1.21e+ 0  2.13e+ 0     0.569  0.600 
##  5 n_teton    <gropd_df [6 × 12]> year   1.90e+ 0  2.58e+ 0     0.736  0.503 
##  6 paintbrush <gropd_df [5 × 12]> year   5.70e+ 0  3.70e+ 0     1.54   0.221 
##  7 s_ak_basin <gropd_df [4 × 12]> year   5.71e- 1  2.05e+ 0     0.279  0.806 
##  8 s_cascade  <gropd_df [3 × 12]> year  -1.09e-15  1.72e-16    -6.35   0.0994
##  9 s_teton    <gropd_df [6 × 12]> year   5.54e- 1  4.53e+ 0     0.122  0.909 
## 10 skillet    <gropd_df [3 × 12]> year   5.79e+ 0  2.85e+ 0     2.03   0.291 
## 11 windcave   <gropd_df [3 × 12]> year   1.55e+ 0  8.25e+ 0     0.188  0.882

Non-significant increases at all sites except south cascade

Visualization:

summer_2.5_withps <- merge(summer_3plus, summer.2.5.lms)

summer.2.5 <- ggplot(summer_2.5_withps, aes(x = year, y = dd_2.5, color = source))+
  geom_point()+
  geom_line()+
  geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
  theme_classic()+
  facet_wrap(~full_name, scales = "free")+
  scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
  geom_text_npc(aes(npcx = 0.9, npcy = 0.1, label = paste( "P = ",round(p.value, 2), sep = "")))+
  labs(x = "Date", y = "No. summer (July-Aug.) days with max temp >2.5 C", color = "Stream Source")+
  theme(legend.position = "bottom", 
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1), 
        axis.title = element_text(size = 11, face = "bold"), 
        axis.text = element_text(size = 8), 
        legend.text = element_text(size = 10), 
        legend.title = element_text(size = 11, face = "bold"), 
        strip.text = element_text(size = 10, face = "bold"))
summer.2.5

## `geom_smooth()` using formula 'y ~ x'

## Warning: Removed 12 rows containing non-finite values (stat_smooth).

## Warning: Removed 12 rows containing missing values (geom_point).

## Warning: Removed 4 row(s) containing missing values (geom_path).

Days over 5 C:

summer.5.lms <- summer_nested %>%
  mutate(model = purrr::map(data, ~lm(dd_5~year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
                     tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
  unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
  filter(term == "year")  #remove intercept term to just get info on year for each site
summer.5.lms

## # A tibble: 11 × 7
## # Groups:   site [11]
##    site       data                term  estimate std.error statistic p.value
##    <chr>      <list>              <chr>    <dbl>     <dbl>     <dbl>   <dbl>
##  1 cloudveil  <gropd_df [3 × 12]> year    2.50       1.44      1.73    0.333
##  2 delta      <gropd_df [4 × 12]> year    0          0       NaN     NaN    
##  3 grizzly    <gropd_df [4 × 12]> year    6.7        5.08      1.32    0.318
##  4 mid_teton  <gropd_df [6 × 12]> year   -0.0286     0.188    -0.152   0.887
##  5 n_teton    <gropd_df [6 × 12]> year    2.09       2.33      0.896   0.421
##  6 paintbrush <gropd_df [5 × 12]> year    9.30       4.66      2.00    0.140
##  7 s_ak_basin <gropd_df [4 × 12]> year    0.357      1.27      0.282   0.805
##  8 s_cascade  <gropd_df [3 × 12]> year   -1.61       2.46     -0.656   0.630
##  9 s_teton    <gropd_df [6 × 12]> year    1.05       4.27      0.245   0.818
## 10 skillet    <gropd_df [3 × 12]> year    6.21       2.35      2.64    0.230
## 11 windcave   <gropd_df [3 × 12]> year    1.68       0.729     2.31    0.260

Non-significant increases at all sites except south cascade, mid teton, and delta (0 all years)

Visualization:

summer_5_withps <- merge(summer_3plus, summer.5.lms)

summer.5 <- ggplot(summer_5_withps, aes(x = year, y = dd_5, color = source))+
  geom_point()+
  geom_line()+
  geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
  theme_classic()+
  facet_wrap(~full_name, scales = "free")+
  scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
  geom_text_npc(aes(npcx = 0.9, npcy = 0.1, label = paste( "P = ",round(p.value, 2), sep = "")))+
  labs(x = "Date", y = "No. summer (July-Aug.) days with max temp >5 C", color = "Stream Source")+
  theme(legend.position = "bottom", 
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1), 
        axis.title = element_text(size = 11, face = "bold"), 
        axis.text = element_text(size = 8), 
        legend.text = element_text(size = 10), 
        legend.title = element_text(size = 11, face = "bold"), 
        strip.text = element_text(size = 10, face = "bold"))
summer.5

## `geom_smooth()` using formula 'y ~ x'

## Warning: Removed 12 rows containing non-finite values (stat_smooth).

## Warning: Removed 12 rows containing missing values (geom_point).

## Warning: Removed 4 row(s) containing missing values (geom_path).

Days over 10 C:

summer.10.lms <- summer_nested %>%
  mutate(model = purrr::map(data, ~lm(dd_10~year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
                     tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
  unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
  filter(term == "year")  #remove intercept term to just get info on year for each site
summer.10.lms

## # A tibble: 11 × 7
## # Groups:   site [11]
##    site       data                term   estimate std.error  statistic p.value
##    <chr>      <list>              <chr>     <dbl>     <dbl>      <dbl>   <dbl>
##  1 cloudveil  <gropd_df [3 × 12]> year   0            0     NaN        NaN    
##  2 delta      <gropd_df [4 × 12]> year   0            0     NaN        NaN    
##  3 grizzly    <gropd_df [4 × 12]> year   6.1 e+ 0     6.01    1.02e+ 0   0.417
##  4 mid_teton  <gropd_df [6 × 12]> year   0            0     NaN        NaN    
##  5 n_teton    <gropd_df [6 × 12]> year   2.01e+ 0     3.20    6.30e- 1   0.563
##  6 paintbrush <gropd_df [5 × 12]> year   5.30e+ 0     3.70    1.43e+ 0   0.247
##  7 s_ak_basin <gropd_df [4 × 12]> year  -5.93e-17     0.327  -1.81e-16   1    
##  8 s_cascade  <gropd_df [3 × 12]> year   0            0     NaN        NaN    
##  9 s_teton    <gropd_df [6 × 12]> year   1.57e+ 0     3.54    4.43e- 1   0.681
## 10 skillet    <gropd_df [3 × 12]> year   0            0     NaN        NaN    
## 11 windcave   <gropd_df [3 × 12]> year   0            0     NaN        NaN

Nothing significant.

Visualization:

summer_10_withps <- merge(summer_3plus, summer.10.lms)

summer.10 <- ggplot(summer_10_withps, aes(x = year, y = dd_10, color = source))+
  geom_point()+
  geom_line()+
  geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
  theme_classic()+
  facet_wrap(~full_name, scales = "free")+
  scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
  geom_text_npc(aes(npcx = 0.9, npcy = 0.1, label = paste( "P = ",round(p.value, 2), sep = "")))+
  labs(x = "Date", y = "No. summer (July-Aug.) days with max temp >10 C", color = "Stream Source")+
  theme(legend.position = "bottom", 
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1), 
        axis.title = element_text(size = 11, face = "bold"), 
        axis.text = element_text(size = 8), 
        legend.text = element_text(size = 10), 
        legend.title = element_text(size = 11, face = "bold"), 
        strip.text = element_text(size = 10, face = "bold"))
summer.10

## `geom_smooth()` using formula 'y ~ x'

## Warning: Removed 12 rows containing non-finite values (stat_smooth).

## Warning: Removed 12 rows containing missing values (geom_point).

## Warning: Removed 4 row(s) containing missing values (geom_path).

TABLE 1: Stream info

Make and save a table with stream names, years of data available for each, and stream source:

First, table set-up:

stream_info <- as.data.frame(summer_3plus) %>%
  dplyr::select(full_name, year, source)%>% #dplyr::select full name, year, and source columns
  group_by(full_name, source)%>% #group by full name and source
  summarise(range = paste(min(year, na.rm = T), max(year, na.rm = T), sep = "-"))%>% #add date range column for each stream
  arrange(source, full_name)%>% #sort by source then stream name
  mutate(source = case_when(source == "glacier"~"Glacier", source == "snowmelt"~"Snow melt", source == "sub_ice"~"Subterranean ice"))

## `summarise()` has grouped output by 'full_name'. You can override using the
## `.groups` argument.

Convert to FlexTable and export:

tbl1 <- flextable(stream_info)%>%
  theme_vanilla()%>%
  align(part = "all", align = "left")%>%
  set_header_labels(full_name = "Stream Name", source = "Hydrologic Source", range = "Date Range")%>%
  set_table_properties(layout = "autofit")
tbl1

Stream Name	Hydrologic Source	Date Range
Cloudveil	Glacier	2019-2021
Delta	Glacier	2017-2021
Middle Teton	Glacier	2015-2021
Skillet	Glacier	2018-2021
Grizzly	Snow melt	2018-2021
North Teton	Snow melt	2015-2021
Paintbrush	Snow melt	2017-2021
South Teton	Snow melt	2017-2021
S. AK Basin	Subterranean ice	2016-2021
South Cascade	Subterranean ice	2015-2021
Wind Cave	Subterranean ice	2015-2020

save_as_docx(tbl1,path = "manuscript/tables/table1_streaminfo.docx", align = "center")
save_as_image(tbl1,path = "manuscript/tables/table1_streaminfo.png", res = 300)

## [1] "manuscript/tables/table1_streaminfo.png"

FIGURE 1: Map of sites used:

#Data set-up:
mean.long <- mean(summer_3plus$long) 
mean.lat <- mean(summer_3plus$lat) #calculate mean lat/long of all sites

map_dat <- as.data.frame(summer_3plus)%>%
  dplyr::select(full_name, source, lat, long)%>%
  distinct()%>%
  filter(!row_number(.)==10)

# Basemap: use get_map to pull map from google API:
basemap <- get_map(c(mean.long,mean.lat), #center map on mean lat and long of all sites
                   source = "stamen", zoom = 11, maptype = "terrain-background") #set map type and zoom

## ℹ <]8;;https://maps.googleapis.com/maps/api/staticmap?center=43.741346,-110.824899&zoom=11&size=640x640&scale=2&maptype=terrain&key=xxxhttps://maps.googleapis.com/maps/api/staticmap?center=43.741346,-110.824899&zoom=11&size=640x640&scale=2&maptype=terrain&key=xxx]8;;>

## ℹ Map tiles by Stamen Design, under CC BY 3.0. Data by OpenStreetMap, under ODbL.

bb <- attr(basemap, "bb") #extract  max/min lat and long of the basemap for clipping

##GTNP boundary:
gtnp <- rgdal::readOGR("gis_data/gnp_gtnp")%>% #read in stored shapefile of GTNP and GNP polygons
  spTransform(crs("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")) #Update CRS to lat/long (comes in as UTM)

## Please note that rgdal will be retired during October 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2023/05/15/evolution4.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-7, (SVN revision 1203)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.4.2, released 2022/03/08
## Path to GDAL shared files: /Users/GordonG/Library/R/x86_64/4.2/library/rgdal/gdal
## GDAL binary built with GEOS: FALSE 
## Loaded PROJ runtime: Rel. 8.2.1, January 1st, 2022, [PJ_VERSION: 821]
## Path to PROJ shared files: /Users/GordonG/Library/R/x86_64/4.2/library/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.6-1
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/GordonG/iCloud Drive (Archive)/Documents/Work/hotaling_data_analyist/tasr/teton_trends/gis_data/gnp_gtnp", layer: "gnp_gtnp"
## with 2 features
## It has 20 fields

## Warning: PROJ support is provided by the sf and terra packages among others

gtnp_clip <- crop(gtnp, extent(bb[1,2], bb[1,4], bb[1,1], bb[1,3]))%>% #Clip GTNP shapefile to the extent of the basemap
  fortify() #fortify to convert from SpatialPolygon to dataframe for plotting

## Regions defined for each Polygons

Mapping:

site.map <- ggmap(basemap)+ #basemap: terrain background (no labels) from google maps platform console, centered on avg lat/long of all sites
  geom_polygon(data = gtnp_clip, aes(x = long, y = lat, group = group), color = "black", linetype = "dashed", fill = "darkgreen", alpha = 0.25, size = 0.6)+ #Adding clipped polygon of GTNP
  geom_point(data = map_dat, aes(x = long, y = lat, fill = source, group = NULL), color = "black", size = 3, pch = 21)+ #add  sites, color code by source 
  labs(color = "Hydrologic Source", fill = "Hydrologic Source")+ #updating legend title
  scale_fill_manual(values = pres.pal, labels = c("Glacier", "Snow melt", "Subterranean ice"))+
  theme_void()+ #removing lat/long from x and y axes 
  theme(legend.title = element_text(size = 9, face = "bold"), #adjusting legend title/text/plot title size, face, and position
        legend.text = element_text(size = 9), 
        legend.position = c(0.84, 0.12), #Moving legend into the plot area
        legend.background = element_rect(fill = "white", color = "black", size = 0.5), #adding box behind legend
        legend.margin = margin(t=5, r=5, b=5, l = 5, unit = "pt") #Increasing margins on legend box
        )+
  ggrepel::geom_label_repel(data=map_dat, aes(x=long, y = lat, label=full_name, fill = source), size =3.25, position = position_dodge(width = 1), fontface = "bold", show.legend = FALSE)+ #adding labels for each stream, color-coded by source
  geom_text_npc(aes(npcx = 0.84, npcy = 0.65, label = "Grand Teton \n National Park"), size = 4.5, fontface = 4, color = "black", hjust = 0.5, alpha = 0.25)+ #add label for Grand Teton NP
  geom_text_npc(aes(npcx = 0.2, npcy = 0.8, label = "Targhee \n National Forest"), size = 4.5, fontface = 4, color = "black", hjust = 0.5, alpha = 0.25) #add lable for Targhee NF
site.map

B. Grouped by source

First, group by source and re-calculate max, mean, min, and degree days:

summer_source <- stream_summer %>%
  group_by(source, year)%>% # group by source type and year
  summarise(t_xbar = mean(summer_xbar), 
            t_max = mean(max_xbar), 
            t_min = mean(min_xbar), #calculate mean daily temp, mean daily min and max
            dd2.5_xbar = mean(dd_2.5), 
            dd5_xbar = mean(dd_5), 
            dd10_xbar = mean(dd_10),#calculate mean no. days with max temp >2.5, 5, and 10 C.
            n_sites = length(unique(site)))%>%  #calculate no sites for each year
  filter(!is.na(year), !is.na(source)
         #, !n_sites < 2
         )%>%
    mutate(source = case_when(source == "glacier"~"Glacier", source == "snowmelt"~"Snow melt", source == "sub_ice"~"Subterranean ice"))

## `summarise()` has grouped output by 'source'. You can override using the
## `.groups` argument.

i. Are there differences in min, max, avg temps between sources?

Visualization:

source.means <- ggplot(summer_source, aes(x = source, y = t_xbar, fill = source))+
  geom_boxplot()+
  theme_classic()+
  scale_fill_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Subterranean Ice"))+
  theme(legend.position = "none")+
  labs(x = "", y = "Mean daily summer (July-Aug) temp., C")
source.max <- ggplot(summer_source, aes(x = source, y = t_max, fill = source))+
  geom_boxplot()+
  theme_classic()+
  scale_fill_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Subterranean Ice"))+
  theme(legend.position = "none")+
  labs(x = "Source", y = "Mean daily max. summer (July-Aug) temp., C")
source.min <- ggplot(summer_source, aes(x = source, y = t_max, fill = source))+
  geom_boxplot()+
  theme_classic()+
  scale_fill_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Subterranean Ice"))+
  theme(legend.position = "none")+
  labs(x = "", y = "Mean daily min. summer (July-Aug) temp., C")
source.means|source.max|source.min

Snow melt streams look to be significantly warmer and have higher max and mins - no surprise there. Check significance of these differences:

max.aov <- aov(t_xbar~source, data = summer_source)
summary(max.aov)

##             Df Sum Sq Mean Sq F value   Pr(>F)    
## source       2  55.74  27.868   22.95 2.73e-05 ***
## Residuals   15  18.21   1.214                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(max.aov)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = t_xbar ~ source, data = summer_source)
## 
## $source
##                                   diff       lwr       upr     p adj
## Snow melt-Glacier           3.68929231  2.036898  5.341687 0.0000981
## Subterranean ice-Glacier   -0.08553797 -1.737933  1.566857 0.9900886
## Subterranean ice-Snow melt -3.77483028 -5.427225 -2.122436 0.0000769

Snowmelt significantly warmer than sub-ice or glacier, no difference between sub-ice and glacier fed.

ii. LMs of mean/min/max over time:

First, visualization - re-organize data for easier plotting:

source_plot <- as.data.frame(summer_source) %>%
  dplyr::select(source, year, t_xbar, t_max, t_min)%>%
  pivot_longer(-c(source, year), names_to = "measure", values_to = "temp")

Plotting:

ylab <- expression(bold("Stream Temperature " (degree*C)))

source.summer <- ggplot(source_plot, aes(x = year, y = temp, color = measure))+
  geom_point()+
  geom_line()+
  geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
  theme_classic()+
  facet_wrap(~source, scales = "free_x")+
  scale_color_manual(values = c("darkred","darkblue", "darkgreen"), labels = c("Max. temp.", "Min. temp.", "Mean temp"))+
  labs(x = "Date", y = ylab, color = "Measurement Type")+
  theme(legend.position = "bottom", 
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1), 
        axis.title = element_text(size = 11, face = "bold"), 
        axis.text = element_text(size = 8), 
        legend.text = element_text(size = 10), 
        legend.title = element_text(size = 11, face = "bold"), 
        strip.text = element_text(size = 10, face = "bold"))
source.summer

## `geom_smooth()` using formula 'y ~ x'

Test for significance:

Nest plot data by source and measurement type:

source_nested <- source_plot%>%
  nest(data = -c(source, measure))

Modeling:

source.lms <- source_nested %>%
  mutate(model = purrr::map(data, ~lm(temp~year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
                     tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
  unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
  filter(term == "year")  #remove intercept term to just get info on year for each site
source.lms

## # A tibble: 9 × 8
##   source           measure data     term   estimate std.error statistic p.value
##   <chr>            <chr>   <list>   <chr>     <dbl>     <dbl>     <dbl>   <dbl>
## 1 Glacier          t_xbar  <tibble> year   0.175       0.0225   7.77    0.00148
## 2 Glacier          t_max   <tibble> year   0.134       0.0514   2.61    0.0594 
## 3 Glacier          t_min   <tibble> year   0.181       0.0283   6.38    0.00309
## 4 Snow melt        t_xbar  <tibble> year   0.135       0.412    0.328   0.760  
## 5 Snow melt        t_max   <tibble> year   0.337       0.554    0.609   0.575  
## 6 Snow melt        t_min   <tibble> year   0.0364      0.306    0.119   0.911  
## 7 Subterranean ice t_xbar  <tibble> year  -0.0517      0.0986  -0.524   0.628  
## 8 Subterranean ice t_max   <tibble> year  -0.000823    0.131   -0.00629 0.995  
## 9 Subterranean ice t_min   <tibble> year  -0.0694      0.107   -0.647   0.553

Significant increases in min, mean, and max summer temp in glacier-fed streams, no change in snowmelt fed streams, non-significant increases in sub_ice streams.

Calculate degrees of warming per decade for each:

warming_amounts <- source.lms %>%
  mutate(deg_decade = 10*estimate)%>%
  dplyr::select(source, measure, deg_decade, p.value)%>%
  rename(Source = source, Measure = measure, Degrees.Warming = deg_decade, P = p.value)%>%
  mutate(Measure = case_when(Measure == "t_xbar"~"Mean", Measure == "t_min"~"Minimum", Measure == "t_max"~"Maximum"))
warming_amounts

## # A tibble: 9 × 4
##   Source           Measure Degrees.Warming       P
##   <chr>            <chr>             <dbl>   <dbl>
## 1 Glacier          Mean            1.75    0.00148
## 2 Glacier          Maximum         1.34    0.0594 
## 3 Glacier          Minimum         1.81    0.00309
## 4 Snow melt        Mean            1.35    0.760  
## 5 Snow melt        Maximum         3.37    0.575  
## 6 Snow melt        Minimum         0.364   0.911  
## 7 Subterranean ice Mean           -0.517   0.628  
## 8 Subterranean ice Maximum        -0.00823 0.995  
## 9 Subterranean ice Minimum        -0.694   0.553

TABLE 2: Deg. warming per decade

Convert to FlexTable and export:

tbl2<-flextable(warming_amounts)%>% #convert to flextable
  theme_vanilla()%>% #update theme to "vanilla"
  set_table_properties(align = "center", layout = "autofit")%>% #center align the table on the screen, autofit columns
  align(part = "all", align = "left")%>%#left align all text
  set_header_labels(Source = "Hydrologic Source", 
                    Measure = "Temperature Type", 
                    Degrees.Warming = "Change per decade, C", 
                    P = "P-value")
tbl2

Hydrologic Source	Temperature Type	Change per decade, C	P-value
Glacier	Mean	1.749728195	0.001477039
Glacier	Maximum	1.340152287	0.059442199
Glacier	Minimum	1.805660683	0.003094522
Snow melt	Mean	1.350701233	0.759600210
Snow melt	Maximum	3.371404238	0.575377043
Snow melt	Minimum	0.364068864	0.911137685
Subterranean ice	Mean	-0.517086607	0.627706676
Subterranean ice	Maximum	-0.008229723	0.995279993
Subterranean ice	Minimum	-0.694147310	0.552638240

save_as_image(tbl2, path = "manuscript/tables/degrees_warming.png", res = 300) #save as .png @ 300dpi

## [1] "manuscript/tables/degrees_warming.png"

save_as_docx(tbl2, path ="manuscript/tables/degrees_warming.docx", align = "center") #save as .png @ 300dpi

Adding these P-values to the plot:

FIGURE 2: Temp trends by source

source_pvals <- merge(source.lms, summer_source)%>% #add model info to source dataframe
  mutate(ypos = case_when(measure == "t_max" ~ 0.95, measure == "t_xbar" ~ 0.9, measure == "t_min" ~ 0.85)) #create column for adjusting P-value label Y position
  


source.summer.pv<- source.summer+
  geom_text_npc(data = source_pvals, aes(npcx = 0.1, npcy = ypos, color = measure, label = paste( "P = ",round(p.value, 3), sep = ""))) # add p value labels to plot
source.summer.pv

## `geom_smooth()` using formula 'y ~ x'

Save:

ggsave("Temperature/plots/source_summer.pdf", source.summer.pv, device = "pdf", units = "in", height = 5, width = 6.5, dpi = "retina")

## `geom_smooth()` using formula 'y ~ x'

iii. LMs of degree days over time:

Re-organize data for easier plotting and modeling:

sourcedd_plot <- summer_source %>%
  dplyr::select(source, year, dd2.5_xbar, dd5_xbar, dd10_xbar)%>%
  pivot_longer(-c(source, year), names_to = "dd_type", values_to = "no_days")

Test for significance:

Nest plot data by source and measurement type:

sourcedd_nested <- sourcedd_plot%>%
  nest(data = -c(source, dd_type))

Modeling:

sourcedd.lms <- sourcedd_nested %>%
  mutate(model = purrr::map(data, ~lm(no_days~year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
                     tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
  unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
  filter(term == "year")  #remove intercept term to just get info on year for each site
sourcedd.lms

## # A tibble: 9 × 8
## # Groups:   source [3]
##   source           dd_type   data     term  estimate std.error statistic p.value
##   <chr>            <chr>     <list>   <chr>    <dbl>     <dbl>     <dbl>   <dbl>
## 1 Glacier          dd2.5_xb… <tibble> year    2.43       1.30      1.86    0.136
## 2 Glacier          dd5_xbar  <tibble> year    1.79       1.30      1.38    0.241
## 3 Glacier          dd10_xbar <tibble> year    0          0       NaN     NaN    
## 4 Snow melt        dd2.5_xb… <tibble> year   -2.95       4.01     -0.735   0.503
## 5 Snow melt        dd5_xbar  <tibble> year   -2.08       3.80     -0.546   0.614
## 6 Snow melt        dd10_xbar <tibble> year   -1.36       3.03     -0.449   0.676
## 7 Subterranean ice dd2.5_xb… <tibble> year   -0.518      3.32     -0.156   0.883
## 8 Subterranean ice dd5_xbar  <tibble> year    1.00       1.04      0.963   0.390
## 9 Subterranean ice dd10_xbar <tibble> year    0.0714     0.169     0.423   0.694

No significant changes in degree days in any group.

Plotting with these P-values:

sourcedd_pvals <- merge(sourcedd.lms, sourcedd_plot)%>% #add model info to source dataframe
  mutate(ypos = case_when(dd_type == "dd2.5_xbar" ~ 0.85, dd_type == "dd5_xbar" ~ 0.9, dd_type == "dd10_xbar" ~ 0.95)) #create column for adjusting P-value label Y position

sourcedd.summer<- ggplot(sourcedd_pvals, aes(x = year, y = no_days, color = dd_type))+
  geom_point()+
  geom_line()+
  geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
  theme_classic()+
  facet_wrap(~source, scales = "free_x")+
  scale_color_manual(values = c("darkred","orange", "darkorange"), labels = c(">10 C", ">2.5 C", ">5 C"))+
  labs(x = "Year", y = "Num. degree days", color = "Degree day threshold")+
  theme(legend.position = "bottom", 
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1), 
        axis.title = element_text(size = 11, face = "bold"), 
        axis.text = element_text(size = 8), 
        legend.text = element_text(size = 10), 
        legend.title = element_text(size = 11, face = "bold"), 
        strip.text = element_text(size = 10, face = "bold"))+
  geom_text_npc(aes(npcx = 0.9, npcy = ypos, color = dd_type, label = paste( "P = ",round(p.value, 3), sep = ""))) # add p value labels to plot
sourcedd.summer

## `geom_smooth()` using formula 'y ~ x'

iv. Lumping glacier & snowmelt:

First, combine glacier and snowmelt streams and recalculate mean, min, and max summer temperatures:

lumped_source <- stream_summer %>%
  mutate(source_lumped = ifelse(source == "sub_ice", "sub_ice", "other"))%>% #add new col with lumped sources
  filter(!is.na(year), !is.na(source))%>% #remove rows with no year or source
  group_by(source_lumped, year)%>% #group by lumped source and year
  summarise(t_xbar = mean(summer_xbar), #calculate mean, min, and max summer temps
            max_xbar = mean(max_xbar), 
            min_xbar = mean(min_xbar))%>%
  pivot_longer(!c(source_lumped, year), names_to = "measure", values_to = "temp") #pivot longer for plotting + analysis

## `summarise()` has grouped output by 'source_lumped'. You can override using the
## `.groups` argument.

Nest data by lumped source and measurement type:

lumped_nested <- lumped_source%>%
  nest(data = -c(source_lumped, measure))

lm for each source + year combo:

lumped.lms <- lumped_nested %>%
  mutate(model = purrr::map(data, ~lm(temp~year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
                     tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
  unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
  filter(term == "year")  #remove intercept term to just get info on year for each site
lumped.lms

## # A tibble: 6 × 8
## # Groups:   source_lumped [2]
##   source_lumped measure  data     term   estimate std.error statistic p.value
##   <chr>         <chr>    <list>   <chr>     <dbl>     <dbl>     <dbl>   <dbl>
## 1 other         t_xbar   <tibble> year   0.220       0.290    0.760     0.490
## 2 other         max_xbar <tibble> year   0.345       0.413    0.835     0.451
## 3 other         min_xbar <tibble> year   0.146       0.204    0.716     0.513
## 4 sub_ice       t_xbar   <tibble> year  -0.0517      0.0986  -0.524     0.628
## 5 sub_ice       max_xbar <tibble> year  -0.000823    0.131   -0.00629   0.995
## 6 sub_ice       min_xbar <tibble> year  -0.0694      0.107   -0.647     0.553

Add p-values to data and plot:

lumped_pvals <- merge(lumped_source, lumped.lms)%>% #add model info to source dataframe
  mutate(ypos = case_when(measure == "t_xbar" ~ 0.9, measure == "max_xbar" ~ 0.95, measure == "min_xbar" ~ 0.85)) #create column for adjusting P-value label Y position

lumped.summer<- ggplot(lumped_pvals, aes(x = year, y = temp, color = measure))+
  geom_point()+
  geom_line()+
  geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
  theme_classic()+
  facet_wrap(~source_lumped, scales = "free_x")+
   scale_color_manual(values = c("darkred","darkblue", "darkgreen"), labels = c("Max. temp.", "Min. temp.", "Mean temp"))+
  labs(x = "Year", y = "Stream temperature, C", color = "Measurement Type")+
  theme(legend.position = "bottom", 
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1), 
        axis.title = element_text(size = 11, face = "bold"), 
        axis.text = element_text(size = 8), 
        legend.text = element_text(size = 10), 
        legend.title = element_text(size = 11, face = "bold"), 
        strip.text = element_text(size = 10, face = "bold"))+
  geom_text_npc(aes(npcx = 0.1, npcy = ypos, color = measure, label = paste( "P = ",round(p.value, 3), sep = ""))) # add p value labels to plot
lumped.summer

## `geom_smooth()` using formula 'y ~ x'

ggsave("Temperature/plots/lumped_summer.pdf", lumped.summer, device = "pdf", units = "in", height = 4, width = 6.5, dpi = "retina")

## `geom_smooth()` using formula 'y ~ x'

2. Modeling stream temperature using SNOTEL data

First, calculate April 1 SWE and summer air temp from snotel data; add to mean temperature dataset:

#Calc means across both sites:
snotel_means <- snotel %>%
  mutate(date1 = as.Date(date, format = "%m/%d/%y"), #recode date column as R date
         airtemp_c = as.numeric(airtemp_c))%>% #recode airtemp as numeric
  group_by(date1)%>% #group by date
  summarise(swe_xbar = mean(swe_cm), #calculate mean SWE and airtemp across both stations
            airtemp_xbar = mean(airtemp_c, na.rm = T))%>%
  mutate(mo = month(date1), yr = year(date1)) #Add month and year cols

## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `airtemp_c = as.numeric(airtemp_c)`.
## Caused by warning:
## ! NAs introduced by coercion

#Calculate max SWE and mean summer airtemp for each year:
snotel_summary <- snotel_means %>% 
  group_by(yr)%>% #group by year
  mutate(swe_max = max(swe_xbar))%>% #extract max SWE and add to new col
  ungroup()%>% #ungroup to get all columns back
  filter(mo == 7 | mo == 8)%>% #extract summer months
  group_by(yr)%>% #group by year
  summarise(airtemp_summer = mean(airtemp_xbar), #calculate mean summer air temp
            swe_max = unique(swe_max))%>% #keep max SWE
  rename(year = yr)

#Merge with long-term temperature averages: 
stream_snotel <- merge(summer_3plus, snotel_summary)

A. Site-by-site

i. Correlations between SWE and Summer Temp:

LMs - first, nest by site:

stream_snotel_nested <- stream_snotel %>%
  nest(data = - site)

lm for each site:

swe.temp.lms <- stream_snotel_nested %>%
  mutate(model = purrr::map(data, ~lm(summer_xbar~swe_max, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
                     tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
  unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
  filter(term == "swe_max")  #remove intercept term to just get info on year for each site
swe.temp.lms

## # A tibble: 11 × 7
##    site       data              term     estimate std.error statistic p.value
##    <chr>      <list>            <chr>       <dbl>     <dbl>     <dbl>   <dbl>
##  1 windcave   <tibble [3 × 14]> swe_max -0.00146    0.0127     -0.115  0.927 
##  2 s_cascade  <tibble [3 × 14]> swe_max -0.0165     0.00937    -1.77   0.328 
##  3 mid_teton  <tibble [6 × 14]> swe_max  0.000784   0.00325     0.241  0.821 
##  4 n_teton    <tibble [6 × 14]> swe_max -0.0561     0.0405     -1.39   0.238 
##  5 s_ak_basin <tibble [4 × 14]> swe_max  0.00414    0.00603     0.687  0.563 
##  6 s_teton    <tibble [6 × 14]> swe_max -0.0522     0.0492     -1.06   0.348 
##  7 delta      <tibble [4 × 14]> swe_max -0.00178    0.00155    -1.15   0.370 
##  8 paintbrush <tibble [5 × 14]> swe_max -0.0634     0.0543     -1.17   0.327 
##  9 skillet    <tibble [3 × 14]> swe_max -0.0305     0.00316    -9.67   0.0656
## 10 grizzly    <tibble [4 × 14]> swe_max -0.0567     0.0840     -0.675  0.569 
## 11 cloudveil  <tibble [3 × 14]> swe_max -0.0166     0.0329     -0.504  0.703

Add P values and merge:

swe_summer_ps <- merge(swe.temp.lms, stream_snotel)

swe.summer.temp <- ggplot(swe_summer_ps, aes(x = swe_max, y = summer_xbar, color = source))+
  geom_point()+
  geom_line()+
  geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
  theme_classic()+
  facet_wrap(~full_name, scales = "free")+
  scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
  geom_text_npc(aes(npcx = 0.9, npcy = 0.9, label = paste( "P = ",round(p.value, 2), sep = "")))+
  labs(x = "Annual maxium SWE, cm", y = "Mean daily July-Aug. stream temp., C", color = "Stream Source")+
  theme(legend.position = "bottom", 
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1), 
        axis.title = element_text(size = 11, face = "bold"), 
        axis.text = element_text(size = 8), 
        legend.text = element_text(size = 10), 
        legend.title = element_text(size = 11, face = "bold"), 
        strip.text = element_text(size = 10, face = "bold"))
swe.summer.temp

## `geom_smooth()` using formula 'y ~ x'

ii. Summer stream temp ~ summer air temp:

lm for each site:

air.temp.lms <- stream_snotel_nested %>%
  mutate(model = purrr::map(data, ~lm(summer_xbar~airtemp_summer, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
                     tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
  unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
  filter(term == "airtemp_summer")  #remove intercept term to just get info on year for each site
air.temp.lms

## # A tibble: 11 × 7
##    site       data              term        estimate std.error statistic p.value
##    <chr>      <list>            <chr>          <dbl>     <dbl>     <dbl>   <dbl>
##  1 windcave   <tibble [3 × 14]> airtemp_su…   0.112     0.273    0.410     0.752
##  2 s_cascade  <tibble [3 × 14]> airtemp_su…  -0.251     0.0964  -2.60      0.233
##  3 mid_teton  <tibble [6 × 14]> airtemp_su…   0.0643    0.0829   0.775     0.481
##  4 n_teton    <tibble [6 × 14]> airtemp_su…  -0.457     1.32    -0.346     0.747
##  5 s_ak_basin <tibble [4 × 14]> airtemp_su…   0.348     0.250    1.39      0.299
##  6 s_teton    <tibble [6 × 14]> airtemp_su…   1.49      4.16     0.359     0.738
##  7 delta      <tibble [4 × 14]> airtemp_su…  -0.214     0.0826  -2.59      0.122
##  8 paintbrush <tibble [5 × 14]> airtemp_su…  -0.0494    5.21    -0.00947   0.993
##  9 skillet    <tibble [3 × 14]> airtemp_su…  -0.346     2.58    -0.134     0.915
## 10 grizzly    <tibble [4 × 14]> airtemp_su…   3.18      5.68     0.559     0.632
## 11 cloudveil  <tibble [3 × 14]> airtemp_su…   1.81      1.04     1.74      0.332

Add P values and plot:

air_summer_ps <- merge(air.temp.lms, stream_snotel)

swe.summer.temp <- ggplot(air_summer_ps, aes(x = airtemp_summer, y = summer_xbar, color = source))+
  geom_point()+
  geom_line()+
  geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
  theme_classic()+
  facet_wrap(~full_name, scales = "free")+
  scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
  geom_text_npc(aes(npcx = 0.9, npcy = 0.9, label = paste( "P = ",round(p.value, 2), sep = "")))+
  labs(x = "Mean daily July-Aug. air temp, C", y = "Mean daily July-Aug. stream temp., C", color = "Stream Source")+
  theme(legend.position = "bottom", 
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1), 
        axis.title = element_text(size = 11, face = "bold"), 
        axis.text = element_text(size = 8), 
        legend.text = element_text(size = 10), 
        legend.title = element_text(size = 11, face = "bold"), 
        strip.text = element_text(size = 10, face = "bold"))
swe.summer.temp

## `geom_smooth()` using formula 'y ~ x'

B. Grouped by source

First, add SNOTEL summary info to mean summer temps grouped by source:

source_snotel <- merge(summer_source, snotel_summary)

i. LMs of mean/min/max temp ~ SWE:

Re-organize data for analysis and plotting:

source_swe <- source_snotel %>%
  dplyr::select(source, swe_max, t_xbar, t_min, t_max)%>% #dplyr::select temperature, source, and swe cols
  pivot_longer(!c(source, swe_max), names_to = "measure", values_to = "temp") #pivot longer for plotting + analysis

Nest data by source and measure:

sourceSwe_nested <- source_swe %>%
  nest(data = -c(source, measure))

LMs for each source and measure:

source.swe.lms <- sourceSwe_nested %>%
  mutate(model = purrr::map(data, ~lm(temp~swe_max, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
                     tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
  unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
  filter(term == "swe_max")  #remove intercept term to just get info on year for each site
source.swe.lms

## # A tibble: 9 × 8
##   source           measure data     term    estimate std.error statistic p.value
##   <chr>            <chr>   <list>   <chr>      <dbl>     <dbl>     <dbl>   <dbl>
## 1 Glacier          t_xbar  <tibble> swe_max -0.00523   0.00895    -0.585  0.590 
## 2 Glacier          t_min   <tibble> swe_max -0.00349   0.00961    -0.363  0.735 
## 3 Glacier          t_max   <tibble> swe_max -0.00689   0.00801    -0.861  0.438 
## 4 Snow melt        t_xbar  <tibble> swe_max -0.0491    0.0355     -1.38   0.239 
## 5 Snow melt        t_min   <tibble> swe_max -0.0417    0.0239     -1.74   0.156 
## 6 Snow melt        t_max   <tibble> swe_max -0.0545    0.0532     -1.02   0.364 
## 7 Subterranean ice t_xbar  <tibble> swe_max  0.0192    0.00656     2.92   0.0431
## 8 Subterranean ice t_min   <tibble> swe_max  0.0231    0.00566     4.07   0.0152
## 9 Subterranean ice t_max   <tibble> swe_max  0.0106    0.0139      0.765  0.487

Add p-values and plot:

sourceSwe_pval<-merge(source_swe, source.swe.lms)%>%
  mutate(ypos = case_when(measure == "t_xbar" ~ 0.9, measure == "t_max" ~ 0.95, measure == "t_min" ~ 0.85)) #create column for adjusting P-value label Y position

source.swe<- ggplot(sourceSwe_pval, aes(x = swe_max, y = temp, color = measure))+
  geom_point()+
  geom_line()+
  geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
  theme_classic()+
  facet_wrap(~source, scales = "free_x")+
   scale_color_manual(values = c("darkred","darkblue", "darkgreen"), labels = c("Max. temp.", "Min. temp.", "Mean temp"))+
  labs(x = "Annual maximum SWE, cm", y = "Stream temperature, C", color = "Measurement Type")+
  theme(legend.position = "bottom", 
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1), 
        axis.title = element_text(size = 11, face = "bold"), 
        axis.text = element_text(size = 8), 
        legend.text = element_text(size = 10), 
        legend.title = element_text(size = 11, face = "bold"), 
        strip.text = element_text(size = 10, face = "bold"))+
  geom_text_npc(aes(npcx = 0.9, npcy = ypos, color = measure, label = paste( "P = ",round(p.value, 3), sep = ""))) # add p value labels to plot
source.swe

## `geom_smooth()` using formula 'y ~ x'

Save:

ggsave("Temperature/plots/source_swe.pdf", source.swe, device = "pdf", units = "in", height = 4, width = 6.5, dpi = "retina")

## `geom_smooth()` using formula 'y ~ x'

ii. LMs of mean/min/max temp ~ Air temp:

Re-organize data for analysis and plotting:

source_air <- source_snotel %>%
  dplyr::select(source, airtemp_summer, t_xbar, t_min, t_max)%>% #dplyr::select temperature, source, and swe cols
  pivot_longer(!c(source, airtemp_summer), names_to = "measure", values_to = "temp") #pivot longer for plotting + analysis

Nest data by source and measure:

sourceAir_nested <- source_air %>%
  nest(data = -c(source, measure))

LMs for each source and measure:

source.air.lms <- sourceAir_nested %>%
  mutate(model = purrr::map(data, ~lm(temp~airtemp_summer, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
                     tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
  unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
  filter(term == "airtemp_summer")  #remove intercept term to just get info on year for each site
source.air.lms

## # A tibble: 9 × 8
##   source           measure data     term    estimate std.error statistic p.value
##   <chr>            <chr>   <list>   <chr>      <dbl>     <dbl>     <dbl>   <dbl>
## 1 Glacier          t_xbar  <tibble> airtem…    0.249     0.221     1.13    0.322
## 2 Glacier          t_min   <tibble> airtem…    0.300     0.219     1.37    0.243
## 3 Glacier          t_max   <tibble> airtem…    0.124     0.229     0.542   0.616
## 4 Snow melt        t_xbar  <tibble> airtem…   -0.690     1.12     -0.616   0.571
## 5 Snow melt        t_min   <tibble> airtem…   -0.731     0.780    -0.937   0.402
## 6 Snow melt        t_max   <tibble> airtem…   -0.437     1.61     -0.272   0.799
## 7 Subterranean ice t_xbar  <tibble> airtem…    0.147     0.315     0.468   0.664
## 8 Subterranean ice t_min   <tibble> airtem…    0.312     0.321     0.972   0.386
## 9 Subterranean ice t_max   <tibble> airtem…   -0.117     0.410    -0.284   0.790

Add p-values and plot:

sourceAir_pval<-merge(source_air, source.air.lms)%>%
  mutate(ypos = case_when(measure == "t_xbar" ~ 0.9, measure == "t_max" ~ 0.95, measure == "t_min" ~ 0.85)) #create column for adjusting P-value label Y position

source.air<- ggplot(sourceAir_pval, aes(x = airtemp_summer, y = temp, color = measure))+
  geom_point()+
  geom_line()+
  geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
  theme_classic()+
  facet_wrap(~source, scales = "free_x")+
   scale_color_manual(values = c("darkred","darkblue", "darkgreen"), labels = c("Max. temp.", "Min. temp.", "Mean temp"))+
  labs(x = "Mean July-Aug. air temp., C", y = "Stream temperature, C", color = "Measurement Type")+
  theme(legend.position = "bottom", 
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1), 
        axis.title = element_text(size = 11, face = "bold"), 
        axis.text = element_text(size = 8), 
        legend.text = element_text(size = 10), 
        legend.title = element_text(size = 11, face = "bold"), 
        strip.text = element_text(size = 10, face = "bold"))+
  geom_text_npc(aes(npcx = 0.1, npcy = ypos, color = measure, label = paste( "P = ",round(p.value, 3), sep = ""))) # add p value labels to plot
source.air

## `geom_smooth()` using formula 'y ~ x'

Save:

ggsave("Temperature/plots/source_air.pdf", source.air, device = "pdf", units = "in", height = 4, width = 6.5, dpi = "retina")

## `geom_smooth()` using formula 'y ~ x'

3. Overall trends in stream temperature over time - Time series analysis

A. Site-by-site

B. Grouped by source

II. Invertebrates:

1. Trends in richness, abundance, and density over time

A. Site-by-site

First, calculate richness, mean density, and total abundance for each site for each year; remove sites with only 1 year of data:

oneYr <- density_source%>%
  group_by(site)%>% #group by site
  summarise(nYrs = length(unique(Year)))%>% #count no. years
  filter(nYrs <3) #extract sites w/ only 1-2 yrs data


invert_means <- density_source %>%
  group_by(stream_code, full_name, site, Year)%>% #group by stream and year
  summarise(richness = length(unique(taxa_lwr)),  #calculate yearly richness, density, and total abundance
            density_xbar = mean(density_xbar), 
            tot_abund = sum(abundance_total), 
            source = unique(source))%>%
  filter(!site%in%oneYr$site) #drop sites with 1-2yrs of data

## `summarise()` has grouped output by 'stream_code', 'full_name', 'site'. You can
## override using the `.groups` argument.

i. LM’s of Richness ~ year at each site:

First, nest data by site:

invert_nested <- invert_means %>%
  nest(data = -site)

LM’s of Richness ~ year for each site:

richness.yr.lm <- invert_nested %>%
  mutate(model = purrr::map(data, ~lm(richness~Year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
                     tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
  unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
  filter(term == "Year")  #remove intercept term to just get info on year for each site
richness.yr.lm

## # A tibble: 11 × 7
## # Groups:   site [11]
##    site       data               term   estimate std.error statistic p.value
##    <chr>      <list>             <chr>     <dbl>     <dbl>     <dbl>   <dbl>
##  1 s_ak_basin <gropd_df [6 × 7]> Year  -1.34e+ 0     0.956 -1.41e+ 0  0.233 
##  2 cloudveil  <gropd_df [3 × 7]> Year   1.50e+ 0     0.866  1.73e+ 0  0.333 
##  3 delta      <gropd_df [5 × 7]> Year  -1.10e+ 0     0.907 -1.21e+ 0  0.312 
##  4 grizzly    <gropd_df [5 × 7]> Year   1.00e- 1     0.790  1.27e- 1  0.907 
##  5 mid_teton  <gropd_df [7 × 7]> Year  -8.93e- 1     0.394 -2.27e+ 0  0.0725
##  6 n_teton    <gropd_df [7 × 7]> Year  -1.54e+ 0     0.886 -1.73e+ 0  0.144 
##  7 paintbrush <gropd_df [5 × 7]> Year  -8.00e- 1     0.542 -1.48e+ 0  0.236 
##  8 s_cascade  <gropd_df [6 × 7]> Year  -2.93e+ 0     0.852 -3.44e+ 0  0.0264
##  9 skillet    <gropd_df [5 × 7]> Year  -1.01e-13     1.25  -8.08e-14  1.00  
## 10 s_teton    <gropd_df [7 × 7]> Year  -1.43e+ 0     0.818 -1.75e+ 0  0.141 
## 11 windcave   <gropd_df [7 × 7]> Year  -7.86e- 1     0.340 -2.31e+ 0  0.0686

Add p-values and plot:

richness_withps <- merge(invert_means, richness.yr.lm)

pres.pal <- c("#C26ED6", "#F89225", "#76D96F")

richness.yr <- ggplot(richness_withps, aes(x = Year, y = richness, color = source))+
  geom_point()+
  geom_line()+
  geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
  theme_classic()+
  facet_wrap(~full_name, scales = "free")+
  scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
  geom_text_npc(aes(npcx = 0.9, npcy = 0.9, label = paste( "P = ",round(p.value, 2), sep = "")))+
  labs(x = "Date", y = "Taxonomic richness", color = "Stream Source")+
  theme(legend.position = "bottom", 
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1), 
        axis.title = element_text(size = 11, face = "bold"), 
        axis.text = element_text(size = 8), 
        legend.text = element_text(size = 10), 
        legend.title = element_text(size = 11, face = "bold"), 
        strip.text = element_text(size = 10, face = "bold"))
richness.yr

## `geom_smooth()` using formula 'y ~ x'

Significant decline in richness at S Cascade, non-significant declines at most other sites

Save:

ggsave("invert_data/plots/richness_year.pdf", richness.yr, device = "pdf", dpi = "retina", units = "in", height = 5, width = 6.5)

## `geom_smooth()` using formula 'y ~ x'

ii. LM’s of Density ~ year at each site:

LM’s of Density ~ year for each site:

dens.yr.lm <- invert_nested %>%
  mutate(model = purrr::map(data, ~lm(density_xbar~Year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
                     tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
  unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
  filter(term == "Year")  #remove intercept term to just get info on year for each site
dens.yr.lm

## # A tibble: 11 × 7
## # Groups:   site [11]
##    site       data               term  estimate std.error statistic p.value
##    <chr>      <list>             <chr>    <dbl>     <dbl>     <dbl>   <dbl>
##  1 s_ak_basin <gropd_df [6 × 7]> Year     12.9      12.2     1.06     0.351
##  2 cloudveil  <gropd_df [3 × 7]> Year     64.6      16.9     3.82     0.163
##  3 delta      <gropd_df [5 × 7]> Year     15.9      41.0     0.389    0.724
##  4 grizzly    <gropd_df [5 × 7]> Year      6.06     44.9     0.135    0.901
##  5 mid_teton  <gropd_df [7 × 7]> Year     42.9      39.5     1.09     0.327
##  6 n_teton    <gropd_df [7 × 7]> Year     -7.10     18.6    -0.382    0.718
##  7 paintbrush <gropd_df [5 × 7]> Year    -12.9     203.     -0.0637   0.953
##  8 s_cascade  <gropd_df [6 × 7]> Year     -2.25     16.2    -0.139    0.896
##  9 skillet    <gropd_df [5 × 7]> Year     32.1      30.6     1.05     0.372
## 10 s_teton    <gropd_df [7 × 7]> Year      3.68      2.05    1.79     0.133
## 11 windcave   <gropd_df [7 × 7]> Year     28.6      19.6     1.46     0.205

Add p-values and plot:

density_withps <- merge(invert_means, dens.yr.lm)

pres.pal <- c("#C26ED6", "#F89225", "#76D96F")

dens.yr <- ggplot(density_withps, aes(x = Year, y = density_xbar, color = source))+
  geom_point()+
  geom_line()+
  geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
  theme_classic()+
  facet_wrap(~full_name, scales = "free")+
  scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
  geom_text_npc(aes(npcx = 0.9, npcy = 0.9, label = paste( "P = ",round(p.value, 2), sep = "")))+
  labs(x = "Date", y = "Mean density, individuals/m2", color = "Stream Source")+
  theme(legend.position = "bottom", 
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1), 
        axis.title = element_text(size = 11, face = "bold"), 
        axis.text = element_text(size = 8), 
        legend.text = element_text(size = 10), 
        legend.title = element_text(size = 11, face = "bold"), 
        strip.text = element_text(size = 10, face = "bold"))
dens.yr

## `geom_smooth()` using formula 'y ~ x'

No significant trends in average density at any site

Save:

ggsave("invert_data/plots/density_year.pdf", dens.yr, device = "pdf", dpi = "retina", units = "in", height = 5, width = 6.5)

## `geom_smooth()` using formula 'y ~ x'

iii. LM’s of total abundance ~ year at each site:

LM’s of Density ~ year for each site:

abund.yr.lm <- invert_nested %>%
  mutate(model = purrr::map(data, ~lm(tot_abund~Year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
                     tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
  unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
  filter(term == "Year")  #remove intercept term to just get info on year for each site
abund.yr.lm

## # A tibble: 11 × 7
## # Groups:   site [11]
##    site       data               term  estimate std.error statistic p.value
##    <chr>      <list>             <chr>    <dbl>     <dbl>     <dbl>   <dbl>
##  1 s_ak_basin <gropd_df [6 × 7]> Year     60.5      149.     0.406   0.705 
##  2 cloudveil  <gropd_df [3 × 7]> Year    333.       133.     2.51    0.242 
##  3 delta      <gropd_df [5 × 7]> Year   -110.       383.    -0.287   0.793 
##  4 grizzly    <gropd_df [5 × 7]> Year    107.       235.     0.455   0.680 
##  5 mid_teton  <gropd_df [7 × 7]> Year      4.89     190.     0.0258  0.980 
##  6 n_teton    <gropd_df [7 × 7]> Year    -28.2      109.    -0.259   0.806 
##  7 paintbrush <gropd_df [5 × 7]> Year   -383.       347.    -1.10    0.350 
##  8 s_cascade  <gropd_df [6 × 7]> Year   -193.        65.0   -2.97    0.0411
##  9 skillet    <gropd_df [5 × 7]> Year    269.       281.     0.955   0.410 
## 10 s_teton    <gropd_df [7 × 7]> Year   -137.        96.4   -1.42    0.215 
## 11 windcave   <gropd_df [7 × 7]> Year     50.6      119.     0.424   0.689

Add p-values and plot:

abund_withps <- merge(invert_means, abund.yr.lm)

pres.pal <- c("#C26ED6", "#F89225", "#76D96F")

abund.yr <- ggplot(density_withps, aes(x = Year, y = tot_abund, color = source))+
  geom_point()+
  geom_line()+
  geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
  theme_classic()+
  facet_wrap(~full_name, scales = "free")+
  scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
  geom_text_npc(aes(npcx = 0.9, npcy = 0.9, label = paste( "P = ",round(p.value, 2), sep = "")))+
  labs(x = "Date", y = "Total abundance", color = "Stream Source")+
  theme(legend.position = "bottom", 
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1), 
        axis.title = element_text(size = 11, face = "bold"), 
        axis.text = element_text(size = 8), 
        legend.text = element_text(size = 10), 
        legend.title = element_text(size = 11, face = "bold"), 
        strip.text = element_text(size = 10, face = "bold"))
abund.yr

## `geom_smooth()` using formula 'y ~ x'

No significant trends in total abundance at any site.

Save:

ggsave("invert_data/plots/abundance_year.pdf", abund.yr, device = "pdf", dpi = "retina", units = "in", height = 5, width = 6.5)

## `geom_smooth()` using formula 'y ~ x'

B. Grouped by source

First, group by source and calculate mean richness, density, and abundance:

source_inverts <- invert_means %>%
  group_by(source, Year)%>%
  summarise(richness_xbar = mean(richness), 
            density_xbar = mean(density_xbar),
            abund_xbar = mean(tot_abund))

## `summarise()` has grouped output by 'source'. You can override using the
## `.groups` argument.

i. LMs of richness ~ year for each source:

First, nest data by source:

source_inverts_nested <- source_inverts %>%
  nest(data = -source)

LMs of richness ~ year for each source:

ricness.source.lm <- source_inverts_nested %>%
  mutate(model = purrr::map(data, ~lm(richness_xbar~Year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
                     tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
  unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
  filter(term == "Year")  #remove intercept term to just get info on year for each site
ricness.source.lm

## # A tibble: 3 × 7
## # Groups:   source [3]
##   source   data             term  estimate std.error statistic p.value
##   <chr>    <list>           <chr>    <dbl>     <dbl>     <dbl>   <dbl>
## 1 glacier  <tibble [7 × 4]> Year    -0.310     0.341    -0.909  0.405 
## 2 snowmelt <tibble [7 × 4]> Year    -1.69      0.895    -1.89   0.117 
## 3 sub_ice  <tibble [7 × 4]> Year    -1.51      0.537    -2.80   0.0379

Add p-values and plot:

richnessSource_pv <- merge(source_inverts, ricness.source.lm)

pres.pal <- c("#C26ED6", "#F89225", "#76D96F")

richness.source <- ggplot(richnessSource_pv, aes(x = Year, y = richness_xbar, color = source))+
  geom_point()+
  geom_line()+
  geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
  theme_classic()+
  facet_wrap(~source)+
  scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
  geom_text_npc(aes(npcx = 0.9, npcy = 0.95, label = paste( "P = ",round(p.value, 2), sep = "")))+
  labs(x = "Year", y = "Mean taxonomic richness", color = "Stream Source")+
  theme(legend.position = "none", 
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1), 
        axis.title = element_text(size = 11, face = "bold"), 
        axis.text = element_text(size = 8), 
        legend.text = element_text(size = 10), 
        legend.title = element_text(size = 11, face = "bold"), 
        strip.text = element_text(size = 10, face = "bold"))
richness.source

## `geom_smooth()` using formula 'y ~ x'

Significant decline in richness at subterranean ice sites

Save:

ggsave("invert_data/plots/richness_source.pdf", richness.source, device = "pdf", dpi = "retina", units = "in", height = 5, width = 6.5)

## `geom_smooth()` using formula 'y ~ x'

ii. LMs of density ~ year for each source:

LMs of density ~ year for each source:

dens.source.lm <- source_inverts_nested %>%
  mutate(model = purrr::map(data, ~lm(density_xbar~Year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
                     tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
  unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
  filter(term == "Year")  #remove intercept term to just get info on year for each site
dens.source.lm

## # A tibble: 3 × 7
## # Groups:   source [3]
##   source   data             term  estimate std.error statistic p.value
##   <chr>    <list>           <chr>    <dbl>     <dbl>     <dbl>   <dbl>
## 1 glacier  <tibble [7 × 4]> Year   -12.1        27.2   -0.444    0.676
## 2 snowmelt <tibble [7 × 4]> Year     0.528      30.4    0.0174   0.987
## 3 sub_ice  <tibble [7 × 4]> Year    10.7        12.2    0.875    0.422

Add p-values and plot:

densSource_pv <- merge(source_inverts, dens.source.lm)

pres.pal <- c("#C26ED6", "#F89225", "#76D96F")

dens.source <- ggplot(densSource_pv, aes(x = Year, y = density_xbar, color = source))+
  geom_point()+
  geom_line()+
  geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
  theme_classic()+
  facet_wrap(~source)+
  scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
  geom_text_npc(aes(npcx = 0.9, npcy = 0.95, label = paste( "P = ",round(p.value, 2), sep = "")))+
  labs(x = "Year", y = "Mean density, individuals/m2", color = "Stream Source")+
  theme(legend.position = "none", 
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1), 
        axis.title = element_text(size = 11, face = "bold"), 
        axis.text = element_text(size = 8), 
        legend.text = element_text(size = 10), 
        legend.title = element_text(size = 11, face = "bold"), 
        strip.text = element_text(size = 10, face = "bold"))
dens.source

## `geom_smooth()` using formula 'y ~ x'

No significant trends in average density in any source.

Save:

ggsave("invert_data/plots/density_source.pdf", dens.source, device = "pdf", dpi = "retina", units = "in", height = 5, width = 6.5)

## `geom_smooth()` using formula 'y ~ x'

iii. LMs of abundance ~ year for each source:

LMs of abundance ~ year for each source:

abund.source.lm <- source_inverts_nested %>%
  mutate(model = purrr::map(data, ~lm(abund_xbar~Year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
                     tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
  unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
  filter(term == "Year")  #remove intercept term to just get info on year for each site
abund.source.lm

## # A tibble: 3 × 7
## # Groups:   source [3]
##   source   data             term  estimate std.error statistic p.value
##   <chr>    <list>           <chr>    <dbl>     <dbl>     <dbl>   <dbl>
## 1 glacier  <tibble [7 × 4]> Year    -103.      159.     -0.649   0.545
## 2 snowmelt <tibble [7 × 4]> Year    -110.       85.0    -1.29    0.254
## 3 sub_ice  <tibble [7 × 4]> Year     -47.4     101.     -0.470   0.658

Add p-values and plot:

abundSource_pv <- merge(source_inverts, abund.source.lm)

pres.pal <- c("#C26ED6", "#F89225", "#76D96F")

abund.source <- ggplot(abundSource_pv, aes(x = Year, y = abund_xbar, color = source))+
  geom_point()+
  geom_line()+
  geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
  theme_classic()+
  facet_wrap(~source)+
  scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
  geom_text_npc(aes(npcx = 0.9, npcy = 0.95, label = paste( "P = ",round(p.value, 2), sep = "")))+
  labs(x = "Year", y = "Mean density, individuals/m2", color = "Stream Source")+
  theme(legend.position = "none", 
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1), 
        axis.title = element_text(size = 11, face = "bold"), 
        axis.text = element_text(size = 8), 
        legend.text = element_text(size = 10), 
        legend.title = element_text(size = 11, face = "bold"), 
        strip.text = element_text(size = 10, face = "bold"))
abund.source

## `geom_smooth()` using formula 'y ~ x'

No significant trends in average abundance in any source.

Save:

ggsave("invert_data/plots/abundance_source.pdf", abund.source, device = "pdf", dpi = "retina", units = "in", height = 5, width = 6.5)

## `geom_smooth()` using formula 'y ~ x'

2. Modeling richness, abundance, and density using SNOTEL data

A. Site-by-site

B. Grouped by source

III. Temperature and Invertebrates

1. Test for correlations between temperature metrics (mean/max/min summer; degree days)

A. Site-by-site

B. Grouped by source